for i=1: NIut
    for j=1:NJut
        if (i<=NIue && j>NJue) continue; end
        
        indexU_cen   = Umap (i  ,j  );
        
        if (j<=NJue)
            
            % Boundary conditions on inlet & outlet
            if (i==1)
                RHS (indexU_cen) = getParabola ( Xp ( indexU_cen,2 ) ) ;
                continue
            end
        else
            % Boundary condition expansion
            if (i==NIue+1)
                RHS (indexU_cen) = 0 ;
                continue
            end
            
        end
        
        if (i==NIut)
            % RHS (indexU_cen) = 0; % Boundary condition oulet - null normal gradient
            U_cen = U (Umap (i,j));
            U_left = U (Umap (i-1,j));
            RHS (indexU_cen) = U_cen - U_cen * (dt/h) * (U_cen-U_left);
            continue
        end
    end
end

for i=1:NIvt
    for j=1:NJvt
        if (i<=NIve && j>NJve) continue; end
        indexV_cen  = dimU + Vmap (i  ,j  ) ;
        
        if (j<=NJve)
            % Boundary conditions on top and bot boundaries
            if (j==1)
                RHS (indexV_cen) = 0;
                continue
            end
            if (j==NJve && i<=NIve )
                RHS (indexV_cen) = 0 ;
                continue
            end
        else
            % Boundary conditions on top and bot boundaries
            if (j==NJvt )
                RHS (indexV_cen) = 0 ;
                continue
            end
        end
        
        
            % Convective outflow modification
        if (i==NIvt)
            U_tr    = U ( Umap ( i+1,j  ) );
            U_tl    = U ( Umap ( i  ,j  ) );
            U_br    = U ( Umap ( i+1,j-1) );
            U_bl    = U ( Umap ( i  ,j-1) );
            V_cen   = V ( Vmap (i,j) ) ;
            V_left  = V ( Vmap ( i-1,j  ) );
            c = (U_tr + U_tl + U_bl + U_br) ; % main convective velocity
            RHS (indexV_cen) = V_cen - (dt*0.25/h) * c * (V_cen-V_left);
            continue
        end
    end
end